Site dilution of quantum spins in the honeycomb lattice 
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We discuss the effect of site dilution on both the magnetization and the density of states of 
quantum spins in the honeycomb lattice, described by the antiferromagnetic Heisenberg spin-S 1 
model. Since the disorder introduced by the dilution process breaks translational invariance, the 
model has to be solved in real space. For this purpose a real-space Bogoliubov-Valatin transformation 
is used. In this work we show that for the S > 1/2 the system can be analyzed in terms of linear 
spin wave theory, in the sense that for all dilution concentrations the assumptions of validity for 
the theory hold. For spin S = 1/2, however, the linear spin wave approximation breaks down. In 
this case, we have studied the effect of dilution on the staggered magnetization using the Stochastic 
Series Expansion Monte Carlo method. Two main results are to be stressed from the Monte Carlo 
method: (i) a better value for the staggered magnetization of the undiluted system, m av (L — > oo) = 
0.2677(6), relatively to the only result available to date in the literature, and based on Trotter 
error extrapolations; (ii) a finite value of the staggered magnetization of the percolating cluster 
at the classical percolation threshold, showing that there is no quantum critical transition driven 
by dilution in the Heisenberg model. In the solution of the problem using linear the spin wave 
method we pay special attention to the presence of zero energy modes. We show, for a finite-size 
system (in a bipartite lattice), that if the two sub-lattices are evenly diluted the system always 
has two zero energy modes, which play the role of Goldstone boson modes for a diluted lattice, 
having no translation symmetry but supporting long range magnetic order. We also discuss the case 
when the two sub-lattices are not evenly diluted. In this case, for finite size lattices, the Goldstone 
modes are not a well defined concept, and special care is needed in taking them into account 
in order for sensible physical results can be obtained. Using a combination of linear spin wave 
analysis and the recursion method we were able to obtain the thermodynamic limit behavior of the 
density of states for both the square and the honeycomb lattices. We have used both the staggered 
magnetization and the density of states to analyze neutron scattering experiments (determining the 
effect of dilution on the system's magnetic moment) and Neel temperature measurements on quasi- 
two-dimensional honeycomb systems. Our results are in quantitative agreement with experimental 
results on Mn p Zni_ p PS3 (a diluted S = 5/2 system) and on the Ba(Ni p Mgi_ p )2V208 (a diluted 
S = 1 system). Our work should stimulate further experimental research in Heisenberg diluted- 
honeycomb systems. 

PACS numbers: 75.10.Jm, 75.50.Ec, 75.30.Ds, 75.40.Mg 



I. INTRODUCTION 

The study of dilution and its effect on the magnetic 
properties of antiferromagnetic materials is a central 
problem in modern condensed matter theoryi 2 ^*^ For 
the square lattice, a number of important experimental 
and theoretical results have been reportedi 3 * 4 *^ For the 
honeycomb lattice, there arc some experimental results in 
the literature^ already, but the corresponding theoretical 
understanding lags far behind. 

Insulating antiferromagnets are possible candidates for 
exhibiting quantum critical points separating ordered 
from disordered phases. The quantum corrections to the 
staggered magnetization of diluted antiferromagnetic in- 
sulators became an important experimental and theoret- 
ical topic after site dilution of La2Cu04 by non mag- 
netic impurities such as Zn or Mg. Theoretical stud- 



ies interpreting the magnetic properties of these diluted 
systems have been recently performedj*^*^ showing a 
good agreement between theory and experiment. A de- 
scription of the effect of dilution on the spin flop phase 
of La2Cu04 was attempted from the point of view of a 
simple mean field theory^ with some qualitative agree- 
ment with experimental results. In addition, the expec- 
tation of a magnetic quantum phase transition driven by 
the interplay of dilution and quantum fluctuations was 
shown not to occur in the antiferromagnetic Heisenberg 
model in a square latticei 4 ^ In the undiluted case, on 
the other hand, it was shown that the Heisenberg model 
itself is incapable of describing the high energy part of 
the spin wave spectrum; a calculation starting from the 
Hubbard model was shown to give the correct high en- 
ergy behavior i^iiSiii 

The key role played by dimensionality in determining 
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the behavior of a system of quantum magnetic moments 
lends special importance to the honeycomb lattice, which 
has the lowest possible co-ordination in more than one 
dimension (see Fig. Realizations of insulating an- 
tiferromagnets based on this lattice have already been 
achieved both with and without magnetic dilution. Re- 
cently Spremo et alH have studied the magnetic proper- 
ties of a metal-organic antiferromagnet on an undiluted 
but distorted honeycomb lattice. The authors found 
good agreement between the theoretical predictions ob- 
tained within the framework of a modified spin wave ap- 
proach and the experimental results for the magnetiza- 
tion as a function of uniform external field and for the 
uniform zero-field susceptibility. 

Honeycomb layers are also found in transition-metal 
thiophosphates MPS3, where M is a first row transition 
metal. These compounds are viewed as "perfect" 2D 
magnetic systems because of the weak van der Waals co- 
hesion energy binding the layers. In each layer the mag- 
netic ions are arranged in a honeycomb lattice. Neutron 
diffraction and magnetic susceptibility studies on M11PS3, 
FePS 3 , and NiPS 3 antiferromagneto 13 i 14 i 15 (5 = 5/2, 
5 = 2, and 5=1, respectively) showed the existence of 
quite different types of ordering among the different com- 
pounds. Whereas for FePS3 and N1PS3 the metal ions are 
coupled ferromagnetically to two of the nearest neigh- 
bors and antiferromagnetically to the third, for M11PS3 
all nearest neighbors interactions within a layer are an- 
tifcrromagnctic. In fact, it turns out that the simplest 
nearest-neighbor antifcrromagnctic Hciscnberg model is 
a reasonable approximation for the description of the 
magnetic properties in MnPS3, although the second- (J2) 
and third- nearest-neighbor (J3) interactions — which are 
both also antiferromagnetic — are not negligible for this 
compound (J1/J2 ~ 10 and J1/J3 ~ 4)>i& Substitu- 
tion of magnetic Mn 2+ ions by nonmagnetic Zn 2+ im- 
purities showed that long-range order (LRO) is lost at 
p = 0.46 ± 0.03 for Mn p Z ni _pPS 3 i 2 i 17 i 18 The fact that 
LRO is preserved for dilutions higher than the classical 
percolation threshold for the honeycomb lattice, p c ~ 0.7, 
is attributed to the significance of J2 and J3 in this com- 
pound. 

Recently, Rogado et alm& have characterized the 
magnetic properties of the 5 = 1 honeycomb com- 
pound BaNi 2 V2 08, which can be described as a weakly 
anisotropic 2D Heisenberg antiferromagnet ^ The mag- 
netic Ni 2+ ions lie on weakly coupled honeycomb lay- 
ers, exhibiting antifcrromagnctic LRO close to 50 K. 
The doped compound Ba(Ni p Mgi_ p )2V208 has a frac- 
tion 1 — p of the honeycomb layer sites substituted by 
Mg 2+ — a nonmagnetic ion. Magnetic susceptibility stud- 
ies showed that the Neel temperature is substantially re- 
duced with increasing doping in the range 0.84 < p < 1. 
For p = 0.84 the onset of antiferromagnetic LRO occurs 
only at Tjy — 17 K, a TV reduction of almost 70% rel- 
ative to its undiluted value. It would be interesting to 
know whether the suppression of antiferromagnetic LRO 
by nonmagnetic impurities occurs at the classical percola- 



tion transition p c ~ 0.7, as predicted by our calculations 
(see below). 

In addition to these exciting experimental results, the 
theoretical result of Mucciolo et al. for the square lattice^ 
where the vanishing of the staggered magnetization for 
the 5 = 1/2 systems coincides with the classical per- 
colation transition, opened the naive expectation that 
for a 2D lattice with nonfrustrating nearest neighbor in- 
teractions and a smaller number of neighbors, magnetic 
quantum-phase transitions driven by the interplay of dis- 
order and quantum fluctuations could occur. The honey- 
comb lattice is the simplest realization of such a lattice, 
for its coordination number is smaller than that of the 
square lattice. On the other hand, large-scale quantum 
Monte Carlo simulations of the square lattice have shown 
that the percolating cluster actually has a robust long- 
range order^ in disagreement with the spin wave calcu- 
lation where this order vanishes very close to the perco- 
lation point. Hence, spin wave theory is not reliable at 
and close to the percolation point for 5 = 1/2, and we 
expect this break-down also for the honeycomb lattice at 
the percolation point. This expectation is confirmed; we 
have performed quantum Monte Carlo simulations that 
show only a rather modest reduction of the sublattice 
magnetization of the percolating cluster, whereas there 
is no long-range order in spin wave theory for 5 = 1/2 in 
this case. 

Experimental 5 = 1/2 antiferromagnetic systems with 
honeycomb lattice structure have already been reported 
by Zhou et alS^ in the A2CuBr4 salt, where A is mor- 
pholinium (C4H10NO). Their data is well described by 
a nearest-neighbor antifcrromagnctic Hciscnberg model, 
but with two different couplings J a and Jb- To the best 
of our knowledge, the dilution of this system has not yet 
been attempted. 

Motivated by the experimental results on diluted 
Mn p Zni_ p PS3 and Ba(NipMgi_p) 2 V208 and by the pos- 
sibility of quantum phase transitions driven by the in- 
terplay of disorder and quantum fluctuations, we study 
here the effect of site dilution on the magnetic proper- 
ties of the Hciscnberg antifcrromagnctic nearest-neighbor 
model, for an arbitrary spin-5 value. Our study is per- 
formed both at zero and finite temperatures. A first at- 
tempt to understand the effect of a nonmagnetic defect 
on the properties of the 5=1/2 2D Heisenberg anti- 
ferromagnet in the honeycomb lattice was made by de 
Chatel et alml In their mean-field approach, a single im- 
purity was introduced in clusters up to 12 spins. It is 
clear, however, that their results can only be applied to 
systems with dilutions up to 1 — p = 1/13. Moreover, the 
random nature of defects cannot be accounted for using 
their method. 

In this paper, we follow the general idea of the work 
of Mucciolo et al.^ by using the linear spin wave ap- 
proximation in real space to compute different physical 
quantities. In addition we use finite-size scaling to deter- 
mine the magnetic moment of the samples. We address 
the problem of determining the density of states (DOS) 
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of our system using a different and more reliable method, 
which gives the behavior of the DOS in the thermody- 
namic limit. The paper is organized as follows: in Sect.lTll 
we present the Hamiltonian and the formalism we use; in 
Sect. Hill we give the numerical details of our method; we 
present the results on the staggered magnetization and 
on density of states as well as on the cluster characteri- 
zation in Sect. Hvl finally, in Sect. we summarize our 
work and present some concluding remarks. 



II. MODEL HAMILTONIAN AND FORMALISM 

The Hciscnberg Hamiltonian describing quantum spins 
in a site-diluted honeycomb lattice is written as 



(i) 



ieA.S 



where S" (S^) is the spin operator on a site i of sublattice 
A (B) . The notation i + S represents a nearest neighbor 
site of site i, connected to i by the vector 8. There are 
three different S vectors given by 



5i = -(l,V3), <5 2 = -(l,-V3) 



-c(l,0), (2) 



where c is the hexagon side length. The rji variables can 
have the values or 1 depending on whether the site i 
exists or not. 

The usual spin wave approximation starts by assuming 
that LRO exists and, in the case of antiferromagnctism, 
that the ground state is not substantially different from 
the Neel state. The mathematical meaning of this simi- 
larity is that the following inequalities should hold: 



5 - (S-' z ) < 5 for i e A, 

5 + (5*' z )<5 forces. 



(3) 
(4) 



With these in mind we express the spin operators in 
terms of bosonic creation and annihilation operators 
as introduced by Holstein and Primakoffi^ Holstcin- 
Primakoff transformation is defined for sublattice A as 




(5) 



In sublattice B the spin have S z = —5 projection in the 
Neel state. Since the bosons should describe excitations 
above the ground state, and this has to be such that in- 
equalities J2J and Q are verified, the S b ' z operator needs 



5 4 b,+ operator must create bosons, and all the operators 
in sublattice B are defined as 

S-' z = -S + b\bi , 



5 



b,- 




(6) 



It is worth mentioning that inequalities and Q can 
also be expressed in terms of the new bosonic operators 
a and b as 



(a\ai) < S for i e A, 
(b\bi) < S for i E B, 



(7) 
(8) 



from which the linear spin wave approximation fol- 
lows straightforwardly by expanding the square roots in 
Eqs. JSJ and ||BJ| in powers of 1/5 and keeping only the 
zcroth order terms: 



5 ? ; — v 25 o>i , 5^ 



25 bi, 
25 b\, 



(9) 



Inserting the resultant approximation 10 into Eq. 
produces the linear spin wave Hamiltonian, which reads 



H = -Jh a S(S + 1) 2J ViVi+s 



i<£A,S 



■JS ^2 rnr]i + s\ha{aia\ + b\ +s bi + s) 
+ a t b l+s + b\ +s a\ 



i£A,S 



(10) 



Note that we have introduced a magnetic anisotropy h a 



m the 5„ a ' 2 5°^ term. 



The linear spin wave Hamiltonian (|10fl can be seen as 
having a classical part of the form 



H cl = -Jh a S(S + 1) 2J ViVi+6 

i£A,8 



(11) 



and a quantum fluctuating part, which can be written as 



Ha, 



({a},{6t}) D ({a},{6t})t, 



(12) 



where ({a},{^})^ is a column vector containing all the 
boson operators and 



D 



K a A 

A T K fc 



(13) 



to be redefined as 5„- 



-5 



b\h. 



is the so-called grand dynamical matrix. For a diluted 
lattice, the number of sites in sublattice A need not be 
the same as that in sublattice B; therefore the dimensions 



Accordingly, the of the blocks in D are N a x N a for K a , Nb X Nb for K 6 
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N a x Nb for A and Nb x N a for A T . The corresponding 
matrix elements arc 



K a 

13 


= JhaSSij-rii^Trii+s , 
s 






for i e A, (14) 


K b - 

13 


= Jh a S5 l jij i ^2 Vi+s , 
s 






for i G 5, (15) 


Aij 


= JSrjiTjj , 


for 


i g A, 


+ (16) 


A ij 


= JSrurij , 


for 


i g B, 


jei + 6. (17) 


The 


diagonalization of 


the 


bosonic Hamiltonian 



amounts to find a transformation T such that 

(T t ) _1 DT~ 1 = diag(^i, . . . ,oj Na ,u> Na+ i, . . .,u Na+Nb ) , 

(18) 

where diag(cJi, . . . , wjv +jv 6 ) stands for a diagonal ma- 
trix with elements wi, . . . , tUN a +N b in its diagonal, N a + 
Nb in number. In this case all the eigenvalues 
u>i, . . . , u>N a , wjVo.+ij • ■ ■ i^Na+Nt are positive. The quasi- 
particles associated with those eigenvalues are obtained 
from 



({a},{/3t})t =T ({a},{&t}) 



t 



(19) 



In the undiluted case, it is well known that Eq. 112(1 can 
be diagonalized through a Bogoliubov-Valatin transfor- 
mation in the reciprocal space. For a subsequent analysis 
it is convenient to reproduce here the results of the cal- 
culation for the undiluted honeycomb latticed We first 
introduce the operators cik and 6k denned as the inverse 
Fourier transforms of dj and , 



1 



'Ok , 



h = 



1 



b k , (20) 



where the k summation ranges over the first Brillouin 
zone of either sublattice A or B. (Do not confuse the 
site index i and the complex imaginary unit also present 
in the Fourier transform). The vector is the position 
vector of site i, and N a = Nb in the absence of dilu- 
tion. Substituting Eq. I120[) into Hamiltonian l|12[) gives 
us H sw = £) k H kl with 



Hv = JSz 



<^kak&-k + 0k & ik a k 



where 0k is defined as 



»ic = - 7 e 
z s 



-ik<5 



, (21) 



(22) 



The diagonalized form of Hamiltonian (|21|l , given by 

#k = ^k(l + a k a k + /? k /?k), (23) 

with 



cdk = JSzy h\ - \4> k \ 



(24) 



can be easily obtained from the following Bogoliubov- 
Valatin transformation. 



«k = M k a k + Vk&Lk > 
/3k = «ka k + "k&-k , 



(25) 



with coefficients itk and «k given as functions of the pa- 
rameters h a and </>k- 

In the diluted case, translational symmetry is lost, and 
the solution in the reciprocal space is as difficult as the 
one in real space. Let us start by constructing an opera- 
tor transformation of the Bogoliubov-Valatin type in real 
space which can be used in the presence of dilution: 



N,, 



OL n = ^2 u m ai + ^2 v nib\ , 

»=1 i=l 
N a N b 



(26) 
(27) 



This definition of cr and @' excitations gives us N a ex- 
type quasi-particles and Nb /3-type quasi-particles. Equa- 
tions 126[) and (|27() define the transformation matrix 11911 , 



U* V* 
W X I ' 



(28) 



where the N a x N a (N b x N a ) and N a x N b (N b x N b ) matri- 
ces U* (W) and V* (X) contain the coefficients u* ni (w n i ) 
and (x n i ) , respectively. Since the quasi-particles must 
have a bosonic character, the quasi-particlc operators 
must satisfy the commutation relations 



„t 1 - 



= [p n ,fl a ]=6 m 



(29) 
(30) 



which lead to the following constraints on the transfor- 
mation coefficients: 



(31) 
(32) 
(33) 



N a 




^ ] u ni u m i 


~ y ] v ni v mi 


i=l 


i=l 




N b 


^2 w ™ w ™ 




i=l 


i=l 




N b 






i=l 


i=l 


N a 


N b 


^ "j w ni u mi 
i=l 


i=l 



(34) 



Equations (|31H34|) can be written in matrix notation 



UU f - VV T = U*U J - V*V J = Iat, 
WW 1 - xxt = W*W T - X*X T 



l N b , 



UW T - VX T = WU T - XV T = , 

w *ut _ x *vt = u*w f - v*x f = , 



(35) 
(36) 
(37) 
(38) 
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where Ijv a (Ii\r 6 ) is the N a x AT a (A^ x Nf,) unit matrix. 
These relations can be put into a more compact form by 
defining the matrix 



-Ijv f 



(39) 



in terms of which Eqs. (|35H38|> can be rewritten as 



Tl p Tt = l p . 



(40) 



Since l p l p = lN a +N b i it can be shown, after simple alge- 
braic transformations, that 

T%T = l p . (41) 

As a result we have four additional (though not indepen- 
dent) sets of orthogonality equations, 



^ ] u ni u n j 
n=l 


- W ni w nj 
n=l 


= h ! 


(42) 


N a 

^ ] v ni v nj 
n=l 


N b 

^ ] x ni x nj 
n=l 




(43) 


N„ 
n=l 


N b 

~ X ni W nj 
n=l 


= o, 


(44) 


N„ 

y ] u ni v nj 
n=l 


N b 

- y ] W ni X n j 
n=l 


= 0. 


(45) 



Since transformations (|26fl and i|27[l diagonalize the spin 
wave Hamiltonian 1|12|) . the quasi-particlcs will obey the 
following commutation relations: 





w n n ' 


(46) 






(47) 




= -<4 a) <4 > 


(48) 






(49) 



[Notice that we have relabeled the positive eigenvalues 
introduced in 118fl to u> [ a ' , . . . , u^' , uJf' , . . . , wj^ .] From 
Eqs. (|46|) and (|26|) we can define an eigenvalue matrix 
equation in the usual form, namely, 



K a —A \ 
A T -K b 



,(<*) 



(50) 



where the column vectors u n and v n contain the coeffi- 
cients u n i and v n i, respectively. From Eq. (|49|) and the 
complex conjugate of Eq. 1|27|) . a similar eigenvalue ma- 
trix equation can be defined, 



K a —A 
A T -K 



Defining the matrices 

fl a = diag(w 1 a) , 



' N a ) i 



(51) 



(52) 



and 



f2^=diag(o;f ) ,...,a;£ ) ) ; 



(53) 



and recalling definition i|28[l for T and definition (|39ll for 
l p , Eqs. H50fl and (|51|l can be combined into a single 
equation, 



(54) 



This can be made more compact still by defining the ma- 
trix VI — diag(wJ Q \ . . . ,cjft ,Wi , ■ ■ ■ ,0Jn ), such that 



Dl„T f = T f m r 



(55) 



From Eq. and the relations l|4Tj|) and i|4"TJl . it is not 
difficult to show that 



lpTlpDlpTUp = CI 



(56) 



Thus, solving the eigenvalue problem defined by Eq. I|54|) 
under the constraint l|40|l is equivalent to finding a trans- 
formation T which satisfies Eq. (|18|) . and where, obvi- 
ously, 



i Th 



(57) 



According to Eq. (|19fl . the diagonalized form of the spin 
wave Hamiltonian is obtained as 

i/ sw ^({a},{&t})D({ a },{6t})t 
= ({a},{/3t})0({a},{/3t})t 

N a N b 



n—1 n—1 

The conclusion of the above discussion is that the op- 
erator transformation given by Eqs. I|26|l and l|27|) di- 
agonalizes the spin wave Hamiltonian (|10|) and that a' 
and /3' are the quasi-particles (with bosonic character) 
associated with the low energy excitations of the antifcr- 
romagnetic Heisenberg Hamiltonian for quantum spins 
in a site diluted honeycomb lattice (needless to say the 
above description applies to other lattices as well). 

Using Eq. (|TT)|) it is possible to write any average of 
the initial bosonic operators in terms of the quasi-particlc 
operators a and (3. The simplest example is the staggered 
magnetization M| tagg at T = given by 



= {N a + N b ) (S - 5m z ) , 



(58) 



where 



N„ 



5m z = J2 Sm{ n ' a *> + Y 6m z l 



(59) 
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with 

^■raS^' (60) 

and 

^-ral^i' (61) 

III. NUMERICAL DETAILS 

The formalism developed in the above section is based 
on the existence of the matrix T and, naturally, on the 
possibility of finding it by some numerical procedure. 
In this work we have used two independent methods 
to compute the transformation matrix and the associ- 
ated eigenenergies. Both methods agree with each other 
within the numerical accuracy of the calculation. One of 
them is based on a Cholesky decomposition and gives the 
T _1 matrix directly, whereas the other solves the eigen- 
value problem defined by Eq. I|54|) , computing the matrix 
TT and from this the matrix T. 



A. Cholesky Decomposition method 

As shown by Colpa^ so long as the grand dynami- 
cal matrix is positive definite, a simple algorithm exists 
for determining T. A hermitian (or symmetric) matrix is 
positive definite if all its eigenvalues are positive. By def- 
inition the quasi-particles and fy have positive or zero 
excitation energy. As will be shown in Subsect. lIII CI the 
zero energy excitations are associated with spin rotations, 
which cost zero energy due to the spin rotational sym- 
metry of the isotropic Heisenberg model. So, provided 
that h a > 1 + , all eigenvalues are positive and the grand 
dynamical matrix is positive definite. The algorithm is 
implemented in three major steps: 

1. for D positive definite a Cholesky decomposition 
can be performed^ and we have D = QQ^, where 
Q is an upper triangular matrix. The existence 
of a Cholesky decomposition guarantees that the 
problem is positive definite; 

2. it can be proved that there exists a unitary 
transformation Y such that Y^Qlp^Q^Y = 
l p ab) diag(wi, . . .,u> Na ,LJ Na+ i, . ..u Na+Nb ); 

3. finally, it can be proved that T _1 = 
Q _1 Y diag(^/ZjT, • ■ • , y/u Na+Nb ). 

B. Bogoliubov-Valatin Transformation method 

The nonhermitian eigcnproblcm defined by Eq. I|54|) 
can be solved with standard numerical algorithms. Here 



we have used subroutines of the LAPACK library. It 
should be noted that the resultant eigenvectors do not 
provide the required matrix directly. After diagonal- 
ization the eigenvectors have to be normalized such that 
they satisfy Eqs. I|31|) and (|32|l for n = m. Degener- 
ate eigenvectors^ have to be carefully analyzed because 
the LAPACK subroutines we have used do not guaran- 
tee that they satisfy Eq. gUJl (though Eqs. to (jSlJ) 
arc satisfied by default for n ^ m). The algorithm is 
implemented as follows: 

1. the matrix Dip' is reduced to an upper Hes- 
senberg form H by an orthogonal transformation 
Y, i.e., H = YDl<, ab) Yt (LAPACK subroutines 
DGEHRD and DORGHR); 

2. the eigenvalues of the upper Hessenberg matrix 
(the same as those of Dip 06 ') and the matrices Q 
and Z from the Schur decomposition H = ZQZ^, 
where Q is an upper quasi-triangular matrix (the 
Schur form), and Z is the orthogonal matrix of 
Schur vectors, are computed (LAPACK subroutine 
DHSEQR); 

3. the right eigenvectors of the upper quasi-triangular 
matrix Q are computed and multiplied by Y^Z^ 
giving the right eigenvectors of T>l p ab ' , whose 
matrix form wc name (LAPACK subroutines 
DTREVC); 

4. degenerate column eigenvectors of are ar- 
ranged in linear combinations such that they satisfy 
Eq. QQ); 

5. nondegenerate column eigenvectors of are nor- 
malized so as to satisfy the orthogonality relations 
of Eqs. (j2U and (|3"2|> . giving matrix T"t; 

6. the positive eigenvalues and respective eigenvectors 
are identified with a modes, and the negative ones 
with the P modes; 

7. eigenvalues and eigenvectors are sorted such that 
matrix has the form defined in Eq. I|28|l . 

C. Zero modes 

It is well known that the clean and isotropic limit of 
Hamiltonian (|10|l has two zero-energy excitations (Gold- 
stone bosons), whose momenta can be determined from 
Eq. H23|) . These gaplcss modes are a consequence of the 
fact that the ground state spontaneously breaks the ro- 
tational symmetry of the Hamiltonian in spin space. It 
can be shownSi that these zero-energy modes have di- 
vergent amplitudes. In two and three dimensions the 
quantum corrections to the staggered magnetization (at 
zero temperature) are finite, meaning that the divergence 
associated with the zero energy modes is intcgrable. We 
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note, however, that if the mean square amplitudes of the 
differences between the two x- and y-components, given 
by 



and 




are computed, we immediately find divergent behavior^ 
The same divergent behavior is also found if we try to get 
the staggered magnetization (|58|l from a finite-size scal- 
ing procedure, doing summations in k-space for finite-size 
systems, including all momenta of the Brillouin zone. 

As shown by Andersonj2L2£ this does not mean that 
the spin wave approximation is breaking down and that 
the system has no LRO. What it means is that these 
divergences are related to the zero point motion of the 
Goldstone modes, and their presence is required to ex- 
ist since in a finite-size system one cannot have solu- 
tions that break the spin rotational symmetry of the 
problem^ The presence of a broken symmetry ground 
state is made possible if we analyze the Hq term in 
Eq. I|21(l . from which the Goldstone modes arise. This 
term cannot be diagonalized through any Bogoliubov- 
Valatin transformation. Actually, it has a continuum 
spectrum starting from the zero energy ground state (see 
Appendix Using this continuum of states we can 
form a wave packet centered around some prefixed orien- 
tation in spin space, with the property of having both a 
finite staggered magnetization, and a mean square roots 
of the quantities (|62(l and (|63|l scaling with 1/N?~ a , with 
a > 0, as long as we pay some extra energy. In addi- 
tion, it can be shown (see Appendix [BJ) that this extra 
energy scales as l/N^ +a , being negligible in the ther- 
modynamic limit. Thus it is a suitable approximation 
to form the above mentioned wave packet from the solu- 
tions of Ho, and to study the energy and the zero point 
motion of all other normal modes within a time interval 
smaller than that needed for the zero-energy wave packet 
to disrupt the coherence of the unidirectional stated The 
understanding of the role played by Goldstone bosons in 
finite-size systems is crucial for computing well defined 
quantities in the thermodynamic limit from calculations 
in finite-size lattices. As a practical example, the stag- 
gered magnetization can be obtained from the finite-size 
calculations if the Goldstone zero point motion contribu- 
tions are subtracted, because the H solutions were al- 
ready used to form the starting broken symmetry state. 
From this procedure we get exactly the same value as 
from the convergent integral in the continuum limit. 

The above discussion now needs to be carried on to 
the diluted case, where the above aspects are more del- 
icate than in the nondisordcred case. In the presence of 



dilution, it is easy to verify that there is at least one zero 
mode in Eq. (|50|) in the isotropic case. This nontrivial 
solution with zero energy satisfies the equation 

m 

with all the amplitudes constant. To prove that this is 
indeed an eigenstate we only need to remember defini- 
tions 1|14H17[) of matrices K a , K b and A, and check that 
the following equalities always hold: 

N a +N b 

Kl = ]T Ay , (65) 

j=N a 

/4=£a£. (66) 

In terms of quasi-particle excitations, the eigenvector 
defined by Eq. I|64|l can be expressed as^ 

N a N b 

al<xJ2 a l+^2 b i> ( 67 ) 

i=l i=l 

in the case of an a-type excitation, and 

N a N b 

AS«E a * + E 6 l< ( 68 ) 
i=i i=i 

if it is a /3-type excitation. Recalling the approximate ex- 
pressions in Eq. (j^J for the operators S < ^ a ' b ' >+ and sl a ' b ^ 
in terms of bosonic operators a and b, Eqs. (|67|l and i|68|) 
can be rewritten as 

a l K s iot > ( 69 ) 
Pi « S+ t . (70) 

Thus, excitations a Q and /?q are precisely the Goldstone 
bosons associated with the broken continuous symmetry 
of spin rotation in the diluted system. 

As will be shown in Subsect. lIII El the thermodynamic 
limit of the staggered magnetization for the diluted sys- 
tem will be obtained from a finite-size scaling analysis. 
As we have started from a broken symmetry ground 
state (the wave packet), which is a direct consequence 
of Eqs. © and (0J, we would proceed as in the clean 
limit and neglect the contributions of ao and Po modes. 
However, although in the undiluted case the number of 
Goldstone modes is always two, when dilution is present 
this number can either be one or two, in a finite size 
lattice. The reason why this is so is that operators S^ t 
and Sfot do not always represent independent excitations, 
i.e., they do not always commute. Naturally S^ ot and S^ ot 
never commute strictly speaking because 

[^tot> ^tot] = 2i?tot ■ (71) 
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Nevertheless, in the clean limit we can easily convince 
ourselves that the expectation value of S* oi is always zero, 
and, as S* ot is a constant of the motion, commutator l|71|l 
will always be zero. To get the value of the commutator 
(|71fl in the presence of dilution we make use of Eq. 
from which one finds 



[^totJ ^tot] 



oc N a - N b . 



(72) 



Now it is easily seen that one can have one or two 
Goldstone modes in a finite-size diluted system: if the 
number of undiluted sites in each sublattice is the same 
(N a = Nb) there will be two Goldstone modes; otherwise, 
if N a ^ Nb, there will be only one. Applying to this 
case the reasoning used for the undiluted case, we should 
then neglect the contributions of the existent Goldstone 
modes. 

As the system size increases, the fluctuations rela- 
tive to the zero mean value of N a — Nb should scale as 
1 / \/N a + Nb, statistically speaking. Therefore the differ- 
ence N a — Nb is again zero in the thermodynamic limit 
and the system has two zero energy excitations. This 
situtation cannot be achieved in finite size lattices, unless 
we restrict ourselvess to cases where the disordered real- 
izations are constrained to obey the condition N a = Nb, 
being clear that the staggered magnetization in the ther- 
modynamic limit cannot depend on this restriction. We 
stress, however, that without this restriction the conclu- 
sions drawn from finite-size lattices would be different 
if we had accepted all sorts of disordered lattice real- 
izations. This difference is due to the contribution of 
the "quasi-divergent" energy modes that emerge when 
N a 7^ Nb. We will get back to this point in Sect. lIVI pre- 
senting numerical evidence for what we have just anal- 
ysed. 



D. Cluster formation and periodic boundary 
conditions 



The study of diluted lattices requires the concept of 
largest cluster, and therefore some care is required in con- 
structing the effective lattice where the quantum problem 
is to be solved. Since we are interested in dilution, the al- 
gorithms discussed in Subsecs. 1111 Al and 1111 Bl arc to be 
implemented not on all occupied lattice sites, but only 
on the sites defined by the largest connected cluster of 
spins, since in the thermodynamic limit a finite magneti- 
zation cannot exist if one is below the percolation critical 
threshold p c . The dilution is induced in the lattice by di- 
luting any site with probability 1 — p. For p = 1 there 
is no dilution at all. When p = p c a classical percolation 
transition occurs in the thermodynamic limit preventing 
the existence of magnetic long range order in the system. 
According to Suding and Ziff^S p c = 0.697043(3) in the 
honeycomb lattice. Here we use p c = 0.697043. 

We work with finite size lattices where periodic bound- 
ary conditions (p.b.c.) are implemented as defined in 
Fig. 2] In Fig. ^ the links on the border are labeled 




FIG. 1: (color online) A finite size honeycomb lattice show- 
ing the periodic boundary conditions used in the numerical 
calculations. 



according to which site they connect to. The lattices 
are characterized by their linear dimension L (L = 3 in 
Fig. nj. The total number of sites for a given L is 2L 2 . 

The algorithm starts by identifying the largest cluster, 
for rigid boundary conditions (this is, with no p.b.c). As 
in Ref . 0, it is only after the largest cluster is found that 
we apply p.b.c. to the original lattice, checking whether 
there are new sites belonging to the largest cluster. As 
previously discussed in Subsect. IIII CI only clusters with 
N a = Nb are to be used, so we reject all disordered lattice 
realizations in which N a ^ Finally, the eigenvalue 

problem is solved for the final cluster using the afore- 
mentioned algorithms. In Fig. [21 we show an example of 
a disorder realization and the corresponding cluster label- 
ing process at p = p c . The larger cluster found for rigid 
boundary conditions can be seen in panel (c) of Fig. 
After p.b.c. implementation the final cluster has a larger 
number of elements (panel (d) in Fig. |5J). 



E. Finite-size scaling 

The eigenvalue problem determines all the eigenvalues 
and eigenfunctions for the cluster, and from these the 
corrections to the staggered magnetization are computed 
according to Eq. (|59fl . For a given p value, N IZ disordered 
lattice realizations with N a = Nb are performed, leading 
to an average staggered magnetization density m av 



1 

~n7-. 



E— — 

i=l 



N % 



(73) 



where Mf ^ [ s the value of Eq. (JSHJ, and N l m is the 
total number of magnetic (undiluted) sites in the lattice, 
for the given disorder realization i. Although m av does 
not depend explicitly on L, the sizes of the clusters are 
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FIG. 2: (color online) An example, for a disorder realization 
in a lattice with L — 10, of the cluster formation. In (a) is 
the original lattice; in (b) each site is chosen to be diluted 
with probability 1 — p c with p c = 0.697043; in (c) the several 
(12 in this case) clusters with rigid boundary conditions have 
been determined; in (d) the larger cluster found in (c) is 
augmented by the periodic boundary conditions. 



determined by L, and therefore different L's lead to dif- 
ferent values for Eq. I|73|l. With this definition we will be 
able to identify m av (p, L — > oo) with the ordered mag- 
netic moment magnitude per magnetic ion measured in 
neutron diffraction experiments. 

From Eq. I|58|) it is easily seen that m av can be ex- 
pressed as the average product of two different contri- 
butions, one purely classical (m l cl ) and the other purely 
quantum (m qm ), 

1 w " 

m av (p, L) = — 2J rn cl m\ m , (74) 
JVrz i=i 

where we used the notation m*j = ~- for the classical 

factor, with TV* = N l a + N l b , and mQ = S - 5m\ for 
the quantum mechanical factor. The quantum contribu- 
tion is simply the staggered magnetization density of the 
larger cluster found in the disorder realization i. It would 
be S in the Neel state but it is reduced by 8m\ (p, L) due 
to quantum fluctuations, whose strength depend on di- 
lution p and lattice size L. If LRO is present we can as- 
sume that the sublatticc magnetization, or equivalcntly 
the staggered magnetization, is a self-averaging quan- 
tity, as was shown to happen in the square lattice case£ 
Thus, in the thermodynamic limit of m a v each disorder 
realization m qm can then be replaced by its infinite-size- 
cxtrapolatcd average, which we denote by m qm , 

m av (p, L — > oo) = m c im qm . (75) 

The classical factor now assumes the standard form for 
the order parameter of the classical percolation problem, 



' (76) 

\ Jv rn / L^oo 

which is zero for p < p c . Therefore a quantum criti- 
cal point can only exist above p c if m qm = for some 
p* > p c - To find m qm we need to compute the aver- 
age infinite-size value of the quantum corrections 5mf 
from our finite size calculations. We show that finite- 
size scaling can be found for this quantity, from which 
results holding in the thermodynamic limit can be ob- 
tained. In our study the size of the largest connected 
cluster N a + Nj, is not fixed, instead the linear dimen- 
sion of the lattice L is. As shown for the square lattice^ 
the alternative approach where the percolating cluster 
size is fixed leads to the same magnetization value in 
the thermodynamic limit. The finite-size scaling proper- 
ties of the quantum correction to the magnetization are 
strictly not known for a disordered system at the perco- 
lation point. However, in practice a direct generalization 
of the pure-system scaling, using the fractal (Hausdorff ) 
dimensionality, has been shown to work well£ Hence we 
will assume 

(8m z ( P , L)) N ^ = 6m? + aL~ D / 2 + bL~ D , (77) 

where Smf is the average quantum correction to the 
staggered magnetization density in the thermodynamic 
limit, and D is the fractal dimension of the cluster, which 
should have the universal value D = 91/48 at p c (in two 
dimensions) , as is confirmed for the square and triangular 
lattices i^i 



F. Density of states 

The real space diagonalization procedure, cither 
Bogoliubov- Valatin or Cholesky decomposition, is very 
time consuming, preventing us from accessing large clus- 
ters (in the honeycomb lattice L = 16 is our upper 
limit). Although for the staggered magnetization den- 
sity a finite-size scaling analysis can be done, we cannot 
easily guess the thermodynamic limit behaviour of the 
DOS from results of systems as small as L = 16. 

In this work the well known recursion method is used to 
compute the average DOS. With this method we can han- 
dle lattices as large as L = 128, with the advantage that 
the obtained DOS is not the typical finite size DOS of a 
system with L = 128, but instead a very good approxi- 
mation for its thermodynamic limit value, guessed from 
this finite-size system. Wc refer the reader to the paper 
of R. Haydock2£ for details in the case of non-interacting 
fermionic systems. Being a real space method the effect 
of disorder can be easily incorporated. Here we adopt 
the formulation introduced in Ref . for disordered elec- 
tronic systems. Further details on the recursion method 
in relation to disordered bosonic bilinear systems [such as 
model Hamiltonian 1)12(1] will be presented elsewhere^ 
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It is worth mentioning that the recursion method has 
proved to be a powerful technique even in the presence 
of interactions^ Actually, the continued fraction repre- 
sentation of the Fourier components of the one particle 
propagator, the basis of the recursion method, is also an 
essential point in the Pade analytical continuation which 
usually arises in the many-body problem^ 

We define the following set of zero temperature re- 
tarded Green's functions in the standard way, 



(*) 



s~iaa 



(t) 



{at(t),b](0)}\0)Q(t), 
{&,(*), aj(0)}|0> 0(f), 

Kt(i), % (o)}|o)e(i), 

{6 i (t),6t(0)}|0)e(t) > 



(78) 



where the notation |0) is used for the ground state of the 
spin wave Hamiltonian (|12(l . The Fourier components of 
each of the Green's functions in Eq. I|78(l . 



G ij (E + iO+) 



are the quantities of interest when determining the DOS. 
Defining the DOS as 



(80) 



it can be easily shown that p(E) is given in terms of the 
Fourier components of the Green's functions l|79|) as 




100 200 300 400 500 600 

iVc 



FIG. 3: (color online) The solid (blue) lines show the distri- 
bution of the number of sites in the largest cluster of a ran- 
domly site-diluted honeycomb lattice. From top to bottom, 
the three panels correspond to dilution levels x = 0.1, 0.5, 0.8. 
From left to right, the peaks correspond to linear system sizes 
L = 4, 5, . . . , 18. The (red) vertical dotted lines indicated the 
average cluster sizes N c , computed as per Eq. 1821 1. 



p{E) = -~Tm 

N r . 7T 



J2Gl a (E + i0 + ) 
]T G^(E + i0+) - ]T G??(-E + i0+) 



iEB 



(81) 



The recursion method gives lm[Gij(E + i0 + )] directly, 
the imaginary part of the Fourier components defined 
in Eq. {75J)ri from which the DOS is straightforwardly 
computed. 



IV. RESULTS 

A. Larger cluster statistics 

The number of sites in a regular planar lattice goes as 
the square of its linear size. In the thermodynamic limit, 
the same scaling applies to the largest cluster of the cor- 
responding randomly-site-diluted lattice. This behaviour 
persists up to the percolation threshold, at which point 
the lattice is dominated by a spanning cluster of frac- 
tal dimension. Beyond percolation, individual clusters 



are no longer extensive: they each constitute a vanishing 
fraction of the total number of sites. 

For a honeycomb lattice of size L and dilution level 
x = (1 — p)/{l — p c ), let P(N c \L,x) denote the probability 
that the largest cluster has N c sites. The average size 
of the largest cluster is simply the corresponding first 
moment: 



2L- 



N c (L,x) = J2 N c P(N c \L,x). 



(82) 



Example probability distributions for the honeycomb lat- 
tice are given in Fig. [3] For small x, the distributions are 
sharply peaked. As a; — > 1, they become progressively 
broader and develop long tails skewed toward small val- 
ues of N c (marking the evolution to a different universal 
scaling function at percolation). 

An effective scaling dimension D e g(L, x) can be de- 
fined by the relation iV c ~ L Dctt . Its evolution with L is 
plotted in Fig. 0] Note that D c g(L, x) has two points of 
attraction in the limit L oo: D c ff(L,x < 1) — ► 2 and 
D c h(L, 1) — > 91/48. Plotted in the appropriate reduced 
coordinates — viz., L D P(N C ) versus L~ D N C where D = 2 
below percolation and D = 91/48 at percolation — the 
probability distribution tends to either a simple delta 
function or the nontrivial curve shown in the inset of 



1 1 




t x=\, D = 91/48 (**) 

' x=l, D = 91/48 

x=l, D = 91/48 (*) 

x=0.9, D = 2 

* = 0.8, D = 2 

x=0.7, D = 2 

- a: = 0.6, D = 2 
> .1 = 0.5, D = 2 

- .t = 0.4, D = 2 
• .1 = 0.3, D = 2 
. .t = 0.2, D = 2 

.t = 0.1, D = 2 

■ .1 = 0, D = 2 

< .1 = 0, RSS 



I 



FIG. 4: (color online) The effective scaling dimension of the 
largest cluster takes one of two values in the L — > oo limit: 
D cB = 2 (0 < x < 1) or D cff = 91/48 (x = 1). For a; < 0.5, 
Z? e fT is close to its asymptotic value at all systems sizes. When 
x is close to 1, very large system sizes are necessary to reach 
the asymptotic regime. The figure inset shows the largest- 
cluster size distribution at percolation plotted in reduced co- 
ordinates. Each curve is computed as a histogram over 10 5 
disorder realizations for system sizes L = 5,6,..., 48. As 
L — > oo, the finite-size results converge to a smooth scaling 
function (one not dissimilar from that of the square-lattice 
case; see Fig. 2 of Ref. 0). 

Fig.H 

As can be seen in Fig.0](inset), a long tail is present for 
smaller cluster sizes. This enhancement of the larger clus- 
ter size distribution can be understood as a consequence 
of the many possible disorder configurations for the same 
dilution. That is, we can have various smaller clusters in- 
stead of one large dominant cluster for the same number 
of diluted sites, though, of course, these disorder config- 
urations are not so favorable. 



B. Finite size scaling analysis 

We have performed numerical real space diagonaliza- 
tion of model Hamiltonian 1|12|) . as described in Sect. 11111 
for the honeycomb and the square lattices. Lattices with 
sizes L = 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 (honeycomb) 
and L = 6,8,10,12,14,16,18,20,22,24,26,28 (square) 
were generated. Averages were taken over N TZ = 10 
disorder realizations^ 

In Figure \5\ we show, for the honeycomb lattice, the 
average quantum correction to the staggered magnetiza- 
tion (Sm z (p, L)) N , for various values of dilution x = (1— 
p)/(l — p c ), as a function of lattice size L~ D ' 2 . The error 
bars are much smaller than the symbols used. The lines 
are fits to the points using the finite-size scaling hypothc- 



FIG. 5: (color online) Finite size scaling of {8m z ) for differ- 
ent values of x — (1— p)/(l — p c ), obtained after 10 5 disorder 
realizations of lattices with equal number of sites in each sub- 
lattice. Also shown for x = 1 is the result obtained when 
the realized lattices are not constrained to have N a = iV;,: 
(*) zero modes were subtracted and the highest amplitude 
(nonzero) mode (see text) was subtracted if N a 7^ Nt; (**) 
only zero modes were subtracted. For x = the RSS result 
was obtained by a reciprocal space sum using the analytical 
result^ 



scs (I77|) . The extrapolated zero abscissa value gives the 
average quantum correction to the staggered magnetiza- 
tion density in the thermodynamic limit 8m J° (p) . In the 
undiluted case there is an excellent agreement between 
the real space diagonalization results (left-triangles) and 
the reciprocal space sum (black squares), obtained from 
the first k-summation in Eq. (C9) of Ref. 0> thus pro- 
viding a reliability test to our algorithms. 

For p = p c we show in Fig. the results obtained 
from three different approaches. The blue up-triangles 
are the results of our standard technique discussed in 
Sect. Hill i.e., only lattices in which A^ a = Nj, were con- 
sidered and zero modes subtracted. The result labeled by 
violet down-triangles refers to a calculation in which the 
disordered realized lattices are not constrained to have 
N a = Nb. The considerable difference between these two 
results is due to the presence of one "quasi-divergent" 
low energy (nonzero) mode when N a ^ N^. That is, even 
though we subtract the zero energy Goldstone mode as 
discussed in Sect. Illll for N a ^ Nb, there is, in this situ- 
ation, a low energy eigenstate that contributes in order 
O(l) for Sm z , compared to the 0(l/N c ) contributions of 
the others eigenstates. If the contribution of this mode is 
subtracted the result labeled by orange diamonds is ob- 
tained, which agrees well with the result of our standard 
technique (where the constrain N a = Nb is always used) . 

To better understand the presence of this nonzero 
energy "quasi-divergent" mode when N a =/= Nb, we 
have computed the contribution to 5m z from the lower 
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FIG. 6: Contributions to Sm z from: the lower nonzero energy 
mode (upper panel); the lower energy mode higher than the 
lower nonzero energy mode (lower panel). The average was 
taken over 10 4 disordered honeycomb lattices, with N a —Nt = 
±1, at p = p c . 

nonzero energy mode (Sm^), and the next one in en- 
ergy (8m), '), constrained to lattices with N a — Nb = ±1. 
Figure H3 shows the behaviour of 5m z (upper panel) 

(2) 

and 5m z (lower panel) with the average cluster size 
N c oc L D . The Sm z 2 ^ contribution decreases with N c , 
signaling the linear increase of the number of modes that 
contribute to 8m z . Instead, the contribution Suiz in- 
creases with N c , and will be of O(l) in the thermody- 
namic limit. As already mentioned in Sect. Hill if N a and 
N b are both of magnitude 10 23 , then, if N a — N b = ±1, 
there will be, for any practical purpose, two Goldstonc 
modes and not only one. This statement should always 
be true if \N a — N b \ <C -/V a ~ Nb. The results presented 
in top panel of Fig.[S]agree with this general picture. Fur- 
thermore, they imply that even for small sizes there is a 
mode, which will be identified with a Goldstone mode in 
the thermodynamic limit, that contributes "macroscopi- 
cally" to 5m z , though having a finite energy. 

C. Staggered magnetization 

The results we found for the quantum mechanical fac- 
tor m qm (x) are summarized in Fig. for the honeycomb 
lattice (panel (a)) and for the square lattice (panel (b)). 
Three different values of spin, S = ^, 1, |, are shown. 

In the undiluted limit we obtain 5m z (0) m 0.258 for 
the honeycomb lattice, and 5m z (0) f» 0.197 for the 
square lattice. These results are in excellent agreement 
with quantum Monte Carlo results, namely, 8m z (0) = 
0.2323(6) for the spin 1/2 Heisenberg antiferromagnet in 
the honeycomb lattice (see Subsect. lIVFfl . and 8m z (0) = 
0.1930(3) in the square latticed 

The effect of the classical factor m c \(x) (not shown) is 



FIG. 7: (color online) Average quantum mechanical factor 
m qm (i) vs dilution x — (1 — p)/(l — p c ) for different values 
of spin the S. Panel (a) shows the results for the honeycomb 
lattice and panel (b) for the square lattice. 



only significant very close to p C7 where it vanishes with 
exponent 5/36iiS Thus, for S > | there is a classically 
driven order disorder transition at p c . For S = \ lin- 
ear spin wave theory predicts a quantum critical point 
in both the honeycomb and square lattices to occur at 
x* = 0.85(1) and x* = 0.98(1), respectively. Similar re- 
sults for the square lattice were obtained in Ref.0, though 
the limited number of averages over disorder prevented 
the authors to distinguish x* from x = 1. 

The predicted quantum critical point is absent in quan- 
tum Monte Carlo calculations, either in the honeycomb 
lattice or in the square latticed As already mentioned in 
Sect, [nj we should not expect the validity of spin wave 
approximation when Sm z ~ S, because inequalities Q 
and (jSJ break down in this situation. This is precisely 
what happens when disorder increases for S = g. 

Comparison with experimental results 

Now we compare our results for the staggered mag- 
netization in the spin wave approximation with avail- 
able experimental measurements on Mn p Zni_ p PS3 and 
Ba(Ni p M gl _ p ) 2 V 2 8 : 

a. MnpZni-pPSs The layered compound MnPS3 is 
a S = 5/2 Heisenberg antiferromagnet This huge spin 
value suggests that the spin wave approximation should 
work well in this case. Indeed, the average magnetic mo- 
ment on the Mn atoms was found to be 4.5(2) \ib at 3.5 
K in the pure material^ in excellent agreement with our 
spin wave result m ss 4.48 (1b- The effect of dilution in 
the average magnetic moment of Mn 2+ ions is presented 
in Fig. |S] Neutron diffraction results on Mn p Zni_ p PS3 
are shown as grey circles^ and the red squares are the 
theoretical results within the linear spin wave approxi- 
mation. To go beyond p c (the first-nearest neighbor per- 
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FIG. 8: (color online) Average magnetic moment per mag- 
netic site as function of dilution p. The linear spin wave result 
for the S — 5/2 Heisenberg antiferromagnet in the honeycomb 
lattice (red squares) is compared with neutron scattering data 
on Mn p Zni_ p PS3 from Ref. (grey circles). 



eolation threshold) we would have to take into account 
second- and third-nearest neighbor couplings in Hamil- 
tonian (JJJ . Nevertheless, the effect of dilution for p < p c 
is already well described by the first-nearest neighbor 
model. Furthermore, the agreement between experimen- 
tal and theoretical results even at p = p c , indicates that 
the primarily effect of second- and third-nearest neigh- 
bor interactions is classical. That is, the existence of one 
largest connected cluster with a finite fraction of spins in 
the thermodynamic limit is guaranteed by this couplings 
for p > p c , but the quantum correction to the staggered 
magnetization density is determined by the smaller first- 
nearest neighbors clusters belonging to this larger one, 
at least for p > p c . Further investigations are needed to 
clarify whether this is the correct pictured 

b. Ba(NipMg\- v )iViO% The layered compound 
BaNi2V208 is a spin 5 = 1 antiferromagnet in a 
honeycomb lattice. Neutron diffraction experiments 
have found, in the pure case, an average magnetic 
moment of 1.55(4) [i B for Ni at 8 KjiS which is in good 
agreement with the spin wave result m sa 1.48 fj. B - 
To our knowledge, the magnetic moment has not yet 
been measured for the diluted compound. Nevertheless, 
the available magnetic susceptibility measurements 
on Ba(Ni p Mgi_p)2V208 for dilutions in the range 
0.84 < p < 1, show that the Neel temperature is strongly 
dependent on the amount of dilutionpiS For the highest 
diluted sample (p = 0.84) a reduction of almost 70% 
relative to the undiluted Neel temperature was found. It 
would be interesting to know whether the suppression of 
antifcrromagnctic LRO by nonmagnetic impurities will 
occur at the classical percolation transition p c ~ 0.7, as 
predicted in our calculations. 



D. Neel Temperature 

The Neel temperature of both Mn p Zni_ p PS3 and 
Ba(Ni p Mgi_p) 2 V208 shows a linear suppression with in- 
creasing dilution 1 — p, 2 *^ a feature that is also seen 
in (quasi-2D) diluted Heisenberg antifcrromagncts with 
square lattice £L2U1 

Within the linear spin wave theory developed in Sees. 
Inland lTlTl for diluted antiferromagnetic systems the finite 
temperature staggered magnetization is given by 



M ^(T) = (j2 S i' Z -JL S i' Z 

\ieA ieB 



(83) 



=N r (S — 8m z — SmT' 



where dm z is the zero-temperature correction to the stag- 
gered magnetization defined in Eq. I|59|l . and SmJ(T) is 
the thermal correction 



(84) 
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with generalized <5mi"' Q ^ and Sm^ 1 '^, 
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and n B (uj) = [e u) l kBT — 1) _1 is the Bose distribution 
function. In the thermodynamic limit the averaged 
over disorder staggered magnetization density can be ex- 
pressed as 



ra av (p, T, L — > oo) = TO c im qm (T), 



(87) 



where m c \ is the classical factor defined in Eq. (|76|l . and 
fnqmiT) is the temperature dependent quantum mechan- 
ical factor, 



x(T) =S-Sr 



,T,L- 



>(T). 



In the undiluted case the thermal correction Sm^^p 
1,T) can be expressed as 



(p=l,T) = 



1 



N a + N b 



E 



h a 



-.n B K), 
(89) 

with w k as in Eq. and k given by Eq. (J55J). The 

summation in k is done in the first Brillouin zone of sub- 
lattice A oy B, and can be replaced by an integration 
when L — ► oo. When h a = 1 the spin wave dispersion 
behaves as oc k in the long wave length limit, sim- 
ilarly to the square lattice case. As a consequence the 
thermal correction to the staggered magnetization de- 
velops a logarithmic divergence, which signals the well 
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known suppression of LRO at T > in the 2D isotropic 
Heisenberg model. 

Therefore, if LRO is present up to Tjy ^ 0, either a 
magnetic anisotropy h a or a finite interplanar exchange 
J± (or both) must be present. If the former is the domi- 
nant effect Tjv can be calculated using the mean-field like 
equation^ 



m qm (T N ) = 0. 



(90) 



In the latter the transition should occur when the inter- 
planar coupling is strong enough to stabilize the LRO in 
comparison with thermal fluctuations: 



J±m 2 {p,T = 0) 



A/2 



k B T N , 



(91) 



The parameter £(p, T) is the inplanc correlation length, 
which characterizes the spin fluctuations of a layered sys- 
tem in a paramagnetic phase. The area of a hexagon 
of side c is given by A = c 2 3\/3/2. The correlation 
length can be calculated in the context of the modified 
spin wave theory^ and in the non-diluted (p = 1) case 
it is exponentially divergent with \/T as T — > 0. The 
mean field picture which leads to Eq. l|9"T]) was proposed 
in Ref. 0, and gives a good description of the variation of 
Tn(p)/Tn(0) with dilution 1 — p in a variety of layered 
compounds with square lattice. 

In the case of MnPS3 a small gap of magnitude AE = 
0.5 meV was found in the spin wave energy at the the 
Brillouin zone center^ This energy gap can be explained 
by either a single-ion anisotropy or a dipole coupling, be- 
ing modeled here by a small magnetic anisotropy h a > 1. 
From the spin wave dispersion l|24[) it is found that 
h a « 1.004 is needed to obtain AE = 0.5 meV (a 
nearest-neighbor exchange of magnitude J = 0.8 meV 
was usedi&). We remark that such a small magnetic 
anisotropy has no effect in the conclusions we have made 
so far based in the isotropic Heisenberg model (h a = 1). 
As an example, the the average magnetic moment on the 
Mn atoms given by spin- wave theory is m ~ 4.48 p,s for 
h a = 1 and m ~ 4.55 /is for h a = 1.004, both in excel- 
lent agreement with the experimental value 4.5(2) at 
3.5 Ki£ Inserting this value of h a into Eq. I|89[) we obtain 
Tjv ~ 70 K as a solution of Eq. H90fl . in agreement with 
the measured value Tjy = 78 K>i4 

Nevertheless a finite interplanar exchange of magni- 
tude Jj_ = 0.0019(2) meV is also present in the MnPS 3 
compound^ With f(p = 1,T N = 78 K) = 27.5 A mea- 
sured by neutron scattering^ and c = 3.5 A^ we ob- 
tain from the mean field equation (|91|l Tm as 6 K. This 
small value of T/v is an indication that the effect of the 
interplanar coupling is not as important as the mag- 
netic anisotropy in stabilizing the LRO. Therefore we use 
Eq. (jnOJ to study the effect of dilution on T N (p). The 
thermal correction 5m^ defined by Eq. (|84|l is computed 
via recursion method (see Subsect. IIII F|) . noting that it 
can be expressed as 



5m T z {T) 



dEn B (E)K{E), 



(92) 



■■ ■■ theory 

♦- ♦ mean field 
O experimental revalues 



0.8 


- ■ U 
m 




■ 


2 0.4 
1- 


S = 5/2 


0.2 





.00 0.90 0.80 0.70 0.60 0.50 

P 



FIG. 9: (color online) rjv(p)/Tjv(0) vs p for S = 5/2. Shown 
are the results obtained by numerically solving Eq. I9UI with 
Srn^ computed applying the recursion method to systems 
with L = 128 and averaging over 200 to 400 disorder realiza- 
tions (squares), the mean-field result of Eq. 1941 (diamonds), 
and experimental results on Mn p Zni_ p PS3 from Ref. (cir- 
cles). 



where the kernel K(E) is given by 



K(E) = — -ilm 

Nr 7T 



+ £ G$(E + 4 0+) + ]T G™(-E + *0+) 

ieB i£A 

+ J2G%(-E + i0+)]. (93) 

It is worth mentioning that with the recursion method 
c5mj can be computed with the same precision (limited 
by the linear size L — 128 of the sample) from the undi- 
luted limit p = 1 to the percolation threshold p = p c . 

The result of numerically solving Eq. I|90|l — with Sm^ 
computed by applying the recursion method to systems 
with L = 128 and averaging over 200 to 400 disorder 
realizations — is shown in Fig. Also shown are the re- 
sults of magnetometry measurements on Mn p Zni_ p PS3 



from Ref. |2|- The difference between the theoretical re- 
sults and experimental values suggests that in opposi- 
tion to the magnetic moment at zero temperature (see 
Fig. |HJ) the effect of second- and third-nearest-neighbor 
couplings should be included to obtain a quantitatively 
correct Neel temperature as dilution is increased. An 
estimation of Tn(p)/Tn(1) can as well be obtained by 
standard mean-field theory, T^ F = lJzS(S + 1)4& Re- 
placing S by the zero temperature staggered magnetiza- 
tion density m av (p) defined in Eq. I|75|) . and assuming 
that the coordination number decreases linearly with di- 
lution, z oc p, the ratio Tn(p)/Tn{1) is given by 



Tn{p) 
Tn(1) 



pm a v{p)[m a v(p) + !]■ 



(94) 
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FIG. 10: (color online) T N (p)/T N (0) vs p for S = 1. Shown 
are the results obtained by numerically solving Eq. 19111 with 
Sm^ computed applying the recursion method to systems 
with L = 128 and averaging over 200 to 400 disorder re- 
alizations (squares), the mean- field result of Eq. |0 (di- 
amonds), and experimental results on Ba(Ni p Mgi_ p )2V208 
from Ref. (circles) . 



In Fig. [§] we show as diamonds the results of Eq. (|94p. 
Although this result reproduces the correct dependence 
on p, it should be stressed that as a mean-field approxi- 
mation the absolute value of T/v(p) is overestimated. 

The effect of dilution on the Neel temperature of 
Ba(NipMgi_j,)2V208 was studied by Rogado et al. for 
dilutions in the range 0.84 < p < The few ex- 

perimental results concerning the magnetic properties of 
BaNi 2 V20s are insufficient to undoubtedly determine the 
model which better describes the magnetic behaviour of 
this compound. Although electron-spin resonance mea- 
surements seem to be well fitted by a weakly anisotropic 
Heisenberg model with easy-plane symmetry (XY), i.e., 
h a < 1 in Hamiltonian UlUfl . the same results can as well 
be explained with the isotropic limit of this models Fur- 
ther experiments would be valuable in determining the 
nature of the LRO observed in this compound, in par- 
ticular inelastic neutron scattering from which the spin 
wave dispersion can be measured. Here we assume that 
a small gap is present at the Brillouin zone center, and 
that it can be modeled by a small uniaxial interaction 
anisotropy, i.e., h a > 1 in Hamiltonian (|1L)|I . In partic- 
ular h a — 1 ~ 10~ 4 is needed to get TV 50 K in the 
undiluted case (a nearest-neighbor exchange of magni- 
tude J»4 meV was used) ~ 

The Tjv(p)/Tjv(1) vs p result obtained by numerically 
solving Eq. (|90|1 for 5 = 1, with 8mF z computed applying 
the recursion method to systems with L = 128 and av- 
eraging over 200 to 400 disorder realizations is shown in 
Fig.llOI(sauares). Also shown are the mean- field result of 
Eq. (|94|> (diamonds) and results of magnetic susceptibil- 
ity data for Ba(Ni p Mg 1 _ p ) 2 V20 8 (circles) (Ref. [ljj. The 
disagreement between the mean- field result (Eq. (1941) ) 



and experimental values can be attributed to the small 
spin S = 1 value, which means higher quantum fluctua- 
tions and less mean-field like behaviour. The theoretical 
result (squares) and the experimental values are in rea- 
sonable agreement, though it seems to worsen as dilution 
increases. It should be noted that the spin-wave theory 
for layered materials is not really adequate at T ~ Tjv, 
and when it is applied to the mean-field like Eq. I|90|) 
it tends to overestimate the absolute value of the Neel 
temperature^! 



E. Density of states 

The effect of dilution has a strong impact on the DOS 
of the system. Since the momentum is no-longer a well 
defined quantum number the spin waves acquire a finite 
lifetime^ As a consequence, the basis that diagonalizcs 
the problem has a very different energy spectrum, which 
implies a different DOS. 

We have calculated the DOS of the antiferromagnetic 
Heisenberg model in the linear spin wave approximation 
for the honeycomb and square lattices in the presence 
of dilution. The recursion method briefly discussed in 
Subsect. 1111 Fl was used to study the variation of the DOS 
with dilution. The method is valid from the undiluted 
p = 1 limit to the percolation threshold p c , and enables 
the access to the whole energy spectrum. The precision 
limit is set by the linear size L of the system, which we 
fix here to L = 128 both in the honeycomb and square 
lattices. 

In Fig. ^2 we show the square lattice DOS at four dif- 
ferent values of dilution x. The depletion of the high 
energy part of the DOS in favor of low energy modes 
is clearly seen as dilution is increased, in agreement 
with the results obtained by exact diagonalizc smaller 
systems^ The two structures visible at around Ej JS = 2 
and 3, which Mucciolo et al~ associated with the break- 
ing of the clean-limit magnon branch into three distinct 
but broad branches, are also evident. 

The DOS for the honeycomb lattice is shown in Fig. 
IT2"1 A decrease in the density of high-frequency states and 
the proportional increase in the density of low-frequency 
ones is also clear as dilution increases. This feature can 
then be viewed as a general effect of the presence of di- 
lution. Structures as those observed in the square lattice 
case, just below E/ JS = 2 and 3, are not so easily iden- 
tified. Nevertheless, a feature of this kind seems to be 
present just below E j JS = 2. To determine whether or 
not it can be associated to the presence of fractons, as in 
the square lattice case^ a more detailed study is needed, 
such as the calculation of the dynamical structure factor 
in the diluted honeycomb lattice. 

The effect of moving spectral weight from the top of 
the band to lower energies due to dilution is accompa- 
nied by the appearance of a set of peaks, starting to 
develop in the high-frequency part of the spectrum for 
small dilution and extending to the entire band as dilu- 
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FIG. 11: (color online) DOS of the Heisenberg antiferromag- 
netic model in the linear spin wave approximation for the 
square lattice. An energy mesh with spacing 0.01 in units of 
JS was used. This results were obtained applying the recur- 
sion method to systems with L — 128, and averaging over 200 
to 400 disorder realizations. The dotted line is the clean limit 
DOS. 



FIG. 12: (color online) DOS of the Heisenberg antiferromag- 
netic model in the linear spin wave approximation for the hon- 
eycomb lattice. An energy mesh with spacing 0.01 in units 
of JS was used. This results were obtained applying the re- 
cursion method to systems with L = 128, and averaging over 
200 to 400 disorder realizations. The dotted line is the clean 
limit DOS. 



tion increases. There is, however, a particular peak that 
deserves special attention. This peak can be seen very 
close to the bottom of the band (E = 0) for x > 0.8 both 
in the honeycomb and square lattice DOS. Figure ^| is a 
zoom of the DOS close to E = at x = x c . Being present 
both in the honeycomb and square lattices, though a bit 
stronger in the former, this peak seems to be a general 
feature associated with dilution. In fact, it is closely re- 
lated to the finitcncss of the quantum corrections to the 
staggered magnetization at zero temperature. 

As shown by Mucciolo et alA, the finiteness of the 
quantum fluctuations reduces to the problem of the con- 
vergence of the integral f* ma " dEp(E)E~ 1 . In Fig. d 
we show a polynomial fit to the low-energy behaviour of 
the DOS (red line in the left side of the peak). Although 
it should be seen as guide to the eyes, we can undoubt- 
edly say that in the low-energy limit the DOS behaves as 
p(E) oc E a with a > 1, and thus the above mentioned in- 
tegral is convergent. This result is consistent with the ex- 
istence of an upper bound for the quantum fluctuations in 
any model with a classically ordered ground state whose 
Hamiltonian can be mapped onto that of a system of 
coupled harmonic oscillators, argued by Mucciolo et al.£ 
This result also agrees with the FSS results presented in 
Subsect. IIIIEI where we found finite values for 8mf(x). 
And the fact that Sm ( ^'(x c ) > 1/2 can be attributed to 
the bad-behaviour of the spin- wave approximation when 
6m z ~ S, as will be shown in the next section. 
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FIG. 13: (color online) Low energy behaviour of the DOS 
of the Heisenberg antiferromagnetic model in the linear spin 
wave approximation for the honeycomb (left) and square 
(right) lattices at x = x c . An energy mesh with spacing 
5 x 10~ 4 in units of JS was used. This results were obtained 
applying the recursion method to systems with L = 128, and 
averaging over 800 disorder realizations. 



F. Quantum Monte Carlo results for 5 = 1/2 

We have performed a Monte Carlo study of the S = 
1/2 quantum Heisenberg antifcrromagnet on the site- 
diluted honeycomb lattice using Stochastic Series Expan- 
sion (SSE)i2L2£i Unlike the spin wave approach described 



17 



in Sects. |H] and IIIII — which should be understood as an 
expansion in the relative reduction of the staggered mo- 
ment 5m z /S — this technique is exact (up to statistical 
uncertainties) and well-behaved even when Sm z ~ S. In 
particular, the SSE Monte Carlo can access the the small- 
s', near-percolation regime where the spin wave calcula- 
tion becomes unreliable. 

We have closely followed the procedure outlined in 
Ref. 0, which treats the site dilution problem on the 
square lattice. To accelerate convergence, we have taken 
advantage of the /3-doubling scheme described therein: 
100 equlibration and 200 sampling sweeps are performed 
at each temperature with the resulting configuration 
(an M-element operator list Sm = [a\, 6J, . . . , [om, &at]) 
used to generate a high-probability initial config- 
uration at the next lowest temperature (S^m = 
[oi, h], . . . , [dM,b M ], [om, bit], ■ ■ ■ , [oi, h]) according to 
the cooling schedule = 2, 4, 8, ... , 2048, 4096. 

A refinement to previous work is that we extrapolate 
the staggered magnetization to the thermodynamic limit 
using two different quantities: 
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FIG. 14: The staggered magnetization of the undiluted 
honeycomb lattice (x = 0, JV C = N = 2L 2 ) is extrapo- 
lated to the thermodynamic limit following Eqs. 196aH and 
lIQfibt . A simultaneous fit of the two data sets yields the value 
m av (L ->(») = 0.2677(6). 



to the corresponding functions on the right-hand side ei- 
ther simulaneously — with parameters m qm , {et^}, {bi} — 



Here, Mf agg = E i£ A - EieB §! * the z-projected f separately-with parameters m qm , {a,} and 



qm i 



staggered magetization and m qm is the quantum mechan- 
ical factor introduced in Subsect. IIII El The notation 
( • )l,x represents an ensemble average over the quantum 
states of the system and over all configurations of the 
size-i lattice with dilution x. The site indices in M| tags 
are understood to range over only the largest connected 
cluster. 

Equation l|95a|) . being linear, is analogous to the quan- 
tity S — Sm z computed via spin wave theory. Equa- 
tion (|95b|) is essentially a structure factor and equivalent 
to Eq. (10) of Ref. H The factors 2 and 3 in Eqs. are 
a consequence of the rotational invariance of the ground 
state. Their particular values follow from the averages 
JdCl \Cl-z\ = 47r/2 and JdCl (Cl ■ zf = An/3 where A is a 
vector ranging over the unit sphere. (Such geometric fac- 
tors are irrelevant to the spin wave case; there the ground 
state is symmetry- broken by explicit construction.) 

As in Ref.H we use the straight-forward generalization 
of the finite-size scaling form for the clean system^ 



{bi}. Verifying that m qm sa m' qm serves as a consistency 
check. 

In the case of the undiluted honeycomb lattice (for 



which m av (L — > oo) 



m 



i), we have simulated lattices 
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[As the discussion in Subsect . II V Al makes clear, this con- 
verges to L~ D I' 2 powerlaw behavior at large L, as in 
Eq. I|77|) .] Numerical measurements of the two quantities 
on the left-hand side of Eqs. (|96a|l and l|96bj> may be fit 



up to linear size L — 32 (i.e., up to 2 x 32 2 = 2048 sites). 
Obscrvables were computed using a bootstrap analysis^ 
of 150 bins of 10 5 samples each (1.5 x 10 6 total Monte 
Carlo sweeps). Best fits to the data, shown in Fig. 1141 
give the thermodynamic limit m av (L — > oo) = 0.2677(6). 
This is somewhat smaller than the square lattice value 
w av (i — > oo) = 0.3070(3)^ a reduction that reflects 
the larger quantum fluctuations on the less meanfield- 
likc honeycomb lattice. 

Note that our value of the staggered magnetization is 
larger than (but consistent with) an earlier Monte Carlo 
measurement due to Reger et al^ (within 1.6 standard 
deviations). It is also, we believe, considerably more ac- 
curate. The Reger group's value of m av (i — > oo) = 
0.22(3) was computed by extrapolating relatively large 
Trotter errors (0.1 < At < 0.2) to At — > and small 
systems sizes (4 < L < 8) to L — > oo. Moveover, their 
analysis supposes that the inverse temperature j3 = 10 is 
sufficiently cold to extract the ground state properties of 
the system, which is very likely incorrect^ 

For the diluted honeycomb lattice, we computed 
the staggered magnetization as an average over 10 5 
randomly-generated disorder realizations. Simulations of 
system sizes up to iV c r* 2000 were extrapolated to the 
thermodynamic limit, as shown in Fig. 1151 The figure 
inset illustrates the dependence of m qm on dilution. 



18 



w 




0.25 



FIG. 15: (color online) The main plot shows an extrapola- 
tion to the thermodynamic limit of twice the z-projected stag- 
gered magnetization for various dilution levels x (as indicated 
by the symbols in the upper-left legend). The lines drawn 
through the data points represent a global fit to Eqs. 196H in 
which m qm (a;), ai(x), bi(x), . . . are treated as powerseries in 
x and varied. The resulting function m qm (x) appears as the 
solid (pink) line in the figure inset alongside Monte Carlo re- 
sults for the square lattice (from Ref. 0) and spinwave results 
for the honeycomb lattice. The (red) errorbars indicate the 
values of m qm extrapolated from each fixed- £ dataset taken 
individually. 



In contrast to the spinwave prediction, we find that 
LRO persists right up to the classical percolation thresh- 
old. The magnitude of the staggered magnetization 
decreases with dilution but does not vanish: m qm = 
0.139(6) at x = 1, which represents a roughly 50% re- 
duction in magnetic moment over the undiluted (x = 0) 
lattice. This is comparable to the effect seen in the square 
lattice where m qm (0) = 0.3070(3) falls to m qm (l) = 
0.150(2). 

We observe that the square- and honeycomb-lattice 
values of m qm are remarkably close in the vicinity of 
x = 1. The likely explanation is that the percolating 
clusters retaining little of the structure of their undi- 
luted parent lattice — are themselves quite similar. Both 
have fractal dimension D = 91/48 and a similar near- 
est neighbour count: with increasing site dilution, the 
average coordination number goes from z hc (0) = 3 and 
g sq (0) = 4 to z hc (l) = 2.22 and z sq (l) = 2.52; see Fig.[H 
The Monte Carlo results arc consistent with our under- 
standing that the quantum fluctuations disrupt the LRO 
in inverse proportion to the number of nearest ncigh- 
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FIG. 16: The disorder-averaged coordination number z(x) is 
plotted as a function of dilution level for infinite square and 
honeycomb lattices. The difference between the two lattice 
types narrows as x — > 1. The undiluted square lattice is 33% 
more coordinated than the honeycomb lattice. At percolation, 
it is only 12% more so. 



bours contributing to the local staggered mean field at 
each site. 



V. SUMMARY AND CONCLUDING REMARKS 

In this work we studied the magnetic properties for 
diluted Hciscnberg models in the honeycomb lattice. Re- 
fined results for the density of states in the square lattice 
case were also reported. We have shown that spin wave 
theory in diluted lattices is quite successful in describ- 
ing the magnetic properties of S > 1/2 systems. On the 
other hand, for S = 1/2, spin wave theory breaks down 
and one has to approach the problem using a Monte Carlo 
method. Contrary to the linear spin wave method the, 
the Monte Carlo method does not allow for the deter- 
mination of the density of states. Having the advantage 
of being rotational invariant by construction, the Monte 
Carlo method does not face the problem of the existence 
of zero energy modes. We have discussed in detail what 
is the physics associated with these modes. In the ther- 
modynamic limit they play the role of Goldstonc modes, 
trying to restore the rotational symmetry of the prob- 
lem, explicitly broken by the spin wave approximation. 
We have shown that in a numerical study these modes 
can not be included in the calculation of operator aver- 
ages, if sensible physical results are to be obtained. This 
is because these modes were already used in the construc- 
tion of the broken symmetry state, as was first discussed 
by P. W. Anderson in his seminal paper on spin waves in 
non-diluted lattices^ 

Our approach allows us to compute both the staggered 
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magnetization and the Neel temperature as function of 
the dilution concentration. In particular, the combina- 
tion of spin wave analysis and the recursion method al- 
lows for the calculation of physical quantities virtually in 
the thermodynamic limit. This possibility was not used 
before in similar studies on the square lattice. 

We have used our results to explain the experi- 
mental data of two Heisenberg honeycomb systems: 
Mn p Zn!„pPS3 (a diluted S = 5/2 system) and 
Ba(Ni p Mg 1 _ p ) 2 V 2 8 (a diluted 5 = 1 system). In 
the first case, the available experimental and theoreti- 
cal studies in the non-diluted regime suggest that second- 
and third-nearest-neighbor interactions play a role on the 
physical properties of the system. This can be seen from 
the fact that the measured magnetic moment of the sam- 
ples is finite beyond the classical site-dilution percolation 
threshold. Our calculation suggests, however, that at low 
temperatures and for p > p c the magnetic moment of 
these samples can be accounted for on the basis of a single 
nearest-neighbor coupling. On the other hand, the cal- 
culation of the Neel temperature using a single nearest- 
neighbor coupling is underestimated, as it should indeed 
be case based on the fact that the magnetic order close 
to the Neel temperature should have a measurable con- 
tribution from the other couplings, which are not much 
smaller than the first nearest-neighbor coupling (the Neel 
temperature for this system using second- and third- 
nearest-neighbor interactions will be studied in a future 
publication). Simple calculations based on simple (Ising 
like) mean filed theories, on the other hand, are very 
much insensitive, by construction, to the microscopic de- 
tails of the system. Therefore, and as long as quantum 
fluctuations are not important, a good agreement with 
the experimental data should be obtained. This is the 
case for Mn p Zni_ p PS3, but not for Ba(Ni p Mgi_ p ) 2 V 2 08 
since its much smaller spin brings about the contribu- 
tions of quantum fluctuations. In the case of the system 
Ba(Ni p Mgi_ p ) 2 V 2 08, there are, unfortunately, no mea- 
surement of its magnetic moment in the diluted phase, 
however, the Neel temperature as function of dilution is 
known from thermodynamic measurements. Our results 
show that in this case, most likely, only the first-nearest- 
neighbor coupling (and a very small magnetic anisotropy) 
are needed to describe the behavior of the Neel temper- 
ature upon dilution. It would be important if further 
investigations on this system could be performed in the 
future. 
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APPENDIX A: DIAGONALIZATION OF # k=0 
The k = term in Hamiltonian 121|) can be expressed 

as 

H = JSz (h a {a a Q + b Q b ) + a b + &o a o) ■ (Al) 

This is a standard bilinear model with two coupled 
modes, which is straightforwardly diagonalized through 
a Bogoliubov-Valatin transformation (Eq. J33J)) when 
h a > 1. In the isotropic h a = 1 case it has an infi- 
nite number of eigenstates with a continuum energy spec- 
trum. 

Let us define the following canonical transformation, 

ao = qi + ipi , (A2) 
b = q 2 + ip2 ■ (A3) 

We use the hat notation to distinguish the operators from 
their eigenvalues. The new generalized "position" q and 
"momentum" p operators satisfy the usual commutation 
relations: 

[a ,4]=l [?i,Pi]=^; (A4) 

[b ,b ]=l => [<72,]3 2 ]=^. (A5) 

After simple algebra we find that the h a = 1 Hamiltonian 
(|A1|) can be written in terms of the new operators <?'s and 
p's as 

H = JSz [(ft + q 2 f + (ft - p 2 ) 2 } . (A6) 

The variables q\ + q 2 and p\ —p 2 can be interpreted as the 
center of mass position and the relative momentum, re- 
spectively, of a two particle system, therefore commuting 
with each other 

[qi + 02, Pi - P-2] = l - - % - = 0. (A7) 

Thus the eigenfunctions of Hamiltonian (|A1(I are given 
in as products of the eigenstates of the operator q\ + q 2 
with eigenstates of the operator p\ — p 2 , 

*q,p(<Zi,9 2 ) = S( qi +q 2 - Q)e^-^\ (A8) 
and the aforementioned continuum spectrum is given by 

E q .p = JSz(Q 2 + P 2 ). (A9) 
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APPENDIX B: ANDERSON BROKEN 
SYMMETRY ANALYSIS 



Thus, to limit (q± — qi) to the range A^_^ 2 we need 



In this appendix we closely follow the ideas developed 
by P. W. Anderson^ to show that the ground state of 
an antiferromagnet should display broken spin rotational 
symmetry, even in the absence of any anisotropy. 

As was shown in Appendix^the operators q\ + qi and 
Pi — f>2 are constants of the motion in a system described 
by the isotropic Hamiltonian (|A1|) . having well defined 
expectation values with zero dispersion. From definitions 
(|A2|) and (|A3fl . and Eqs. J5UJ) and ©, it can be easily 
seen that 



1 



qi + q2 = 



QJ1 



pi - p2 = vAm slt > 



(Bl) 
(B2) 



which explains the constant of motion character of the 
<?i + 92 and pi — pi operators (Stot is a constant of mo- 
tion of the original Hcisenberg model). The uncertainty 
relation ensures us that their canonical conjugates coun- 
terparts will have divergent dispersion. As for Eq. i|BT1) 
and (|B2|) it is not difficult to show that the canonical 
conjugates of p\ — pi and q\ + q\ are, respectively, 



<7i - 92 



V^SNa 



iEB 



b,y 



(B3) 
(B4) 



We want to know how much energy is needed to form 
a wave packet with states <|A8|1 (above the ground state), 
such as the expectation values of operators q\ — qi and 
Pi + f>2 have finite dispersion. Let us limit the fluctu- 
ations of the expectation value (qi — qi) to the range 
Aqj-^. From the uncertainty relation the expectation 
value of p\ — pi must now have a nonzero dispersion, 
whose magnitude is given by 



_L 



(B5) 



lira 



JSz 
4A? . 

91-92 



(B6) 



relatively to the ground state energy. Analogously, to 
limit {pi + pi) to the range A^+p.-, we need 



JSz 
4 A? ~ 

P1+P2 



(B7) 



Defining the the mean square amplitudes of the quanti- 
ties and 



{2N, 



we find from (IB3I and (IB4H that 



(B8) 
(B9) 



a 2 = 2JV» 2 

91— 92 o x ' 

A 2 = 2Ngn 

P1+P2 g y ' 



(BIO) 
(BH) 



(note that q\ — q\ and p\ + j>2 have zero expectation 
value). Inserting (|B10|) and (|B11|) in Eqs. l|B6|) and l|B7|) 

it can be seen that to limit a x or a y to a finite value 
we only need an excess energy of magnitude 1 /N a , and 
hence negligible in the thermodynamic limit. As pointed 
out by Anderson, we can even limit a x or <j y to values 

of magnitude l/Na +a , with a > 0, requiring no energy 
when N a — > oo. 



1 O. P. Vajk, P. K. Mang, M. Greven, P. M. Gehring, and 
J. W. Lynn, Science 295, 1691 (2002). 

2 D. J. Goossens, A. J. Studer, S. J. Kennedy, and T. J. 
Hicks, J. Phys.: Condens. Matter 12, 4233 (2000). 

3 A. L. Chernyshev, Y. C. Chen, and A. H. C. Neto, Phys. 
Rev. Lett. 87, 067209 (2001). 

4 E. R. Mucciolo, A. H. C. Neto, and C. Chamon, Phys. Rev. 
B 69, 214424 (2004). 

5 A. W. Sandvik, Phys. Rev. B 66, 024418 (2002). 

6 S. W. Cheong, A. S. Cooper, L. W. Rupp, B. Batlogg, J. D. 
Thompson, and Z. Fisk, Phys. Rev. B 44, 9739 (1991). 

7 A. L. Chernyshev, Y. C. Chen, and A. H. C. Neto, Phys. 



Rev. B 65, 104407 (2002). 

8 N. M. R. Peres, J. Phys.: Condens. Matter 15, 7271 
(2003). 

9 N. M. R. Peres and M. A. N. Araiijo, Phys. Stat. Sol. B 
236, 523 (2003). 

N. M. R. Peres and M. A. N. Araiijo, Phys. Rev. B 65 
(2002). 

1 P. Sengupta, R. T. Scalettar, and R. R. P. Singh, Phys. 
Rev. B 66 (2002). 

2 I. Spremo, F. Schuetz, P. Kopietz, V. Pashchenko, B. Wolf, 
M. Lang, J. W. Bats, C. Hu, and M. U. Schmidt, cond- 
mat/0505425. 



21 



13 G. L. Flem, R. Brec, G. Ouvard, A. Louisy, and P. Seg- 
ransan, J. Phys. Chem. Solids 43, 455 (1982). 

14 K. Kurosawa, S. Saito, and Y. Yamaguchi, J. Phys. Soc. 
Jpn. 52, 3919 (1983). 

15 P. A. Joy and S. Vasudevan, Phys. Rev. B 46, 5425 (1992). 

16 A. R. Wildes, B. Roessli, B. Lebech, and K. W. Godfrey, 
J. Phys.: Condens. Matter 10, 6417 (1998). 

17 N. Chandrasekharan and S. Vasudevan, Phys. Rev. B 54, 
14903 (1996). 

18 D. J. Goossens and T. J. Hicks, J. Phys.: Condens. Matter 
10, 7643 (1998). 

19 N. Rogado, Q. Huang, J. W. Lynn, A. P. Ramirez, D. Huse, 
and R. J. Cava, Phys. Rev. B 65, 144443 (2002). 

20 M. Heinrich, H. K. Von Nidda, A. Loidl, N. Rogado, and 
R. J. Cava, Phys. Rev. Lett. 91 (2003). 

21 P. Zhou and J. E. Drumheller, J. Appl. Phys. 69, 5804 
(1991). 

22 P. F. De Chatel, J. Chadwick, A. M. Mulders, and T. J. 
Hicks, Physica B 344, 117 (2004). 

23 T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940). 

24 N. M. R. Peres, M. A. N. Araiijo, and D. Bozi, Phys. Rev. 
B 70, 195122 (2004). 

25 J. H. P. Colpa, Physica A 93, 327 (1978). 

26 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. 
Flannery, Numerical Recipes in Fortran 77 (Cambridge 
University Press, Cambridge, 1992), p. 89, 2nd ed. 

27 P. W. Anderson, Concepts in Solids (W. A. Benjamin, 
Massachusetts, 1963). 

28 P. W. Anderson, Phys. Rev. 85, 714 (1952). 

29 L. Hulthen, Arkiv Mat. Astron. Fysik 26, 1 (1938). 

30 P. N. Suding and R. M. Ziff, Phys. Rev. E 60, 275 (1999). 

31 L. Pietronero and A. Stella, Physica A 170, 64 (1990). 

32 R. Haydock, in Solid State Physics, edited by H. Ehrenre- 
ich, F. Seitz, and D. Turnbull (Academic Press, New York, 
1980), vol. 35, p. 215. 

33 R. Haerle, R. Haydock, and R. L. Te, Comput. Phys. Com- 
mun. 90, 81 (1995). 

34 E. V. Castro, (to be published). 

35 R. Haydock, Phys. Rev. B 61, 7953 (2000). 

36 K. S. D. Beach, R. J. Gooding, and F. Marsiglio, Phys. 
Rev. B 61, 5147 (2000). 

37 A. W. Sandvik, Phys. Rev. B 56, 11678 (1997). 

38 D. Stauffer and A. Aharony, Introduction to Percolation 
Theory (Taylor and Francis, London, 1992). 

39 E. V. Castro and N. M. R. Peres, (in preparation). 

40 P. Carretta, A. Rigamonti, and R. Sala, Phys. Rev. B 55, 
3734 (1997). 

41 M. Hucker, V. Kataev, J. Pommer, J. Harrass, A. Hosni, 
C. Pfhtsch, R. Gross, and B. Buchner, Phys. Rev. B 59, 
R725 (1999). 

42 K. Takeda, O. Fujita, M. Hitaka, M. Mito, T. Kawae, 
Y. Higuchi, H. Deguchi, Y. Muraoka, K. Zenmyo, H. Kubo, 



et al., J. Phys. Soc. Jpn. 69, 3696 (2000). 

43 M. Takahashi, Phys. Rev. B 40, 2494 (1989). 

44 H. M. Ronnow, A. R. Wildes, and S. T. Bramwell, Physica 
B 276, 676 (2000). 

45 K. Okuda, K. Kurosawa, S. Saito, M. Honda, Z. Yu, and 
M. Date, J. Phys. Soc. Jpn. 55, 4456 (1986). 

46 W. Jones and N. H. March, Theoretical Solid State Physics 
(Dover, New York, 1985), vol. 1, p. 367. 

47 V. Y. Irkhin, A. A. Katanin, and M. I. Katsnelson, Phys. 
Rev. B 60, 1082 (1999). 

48 A. W. Sandvik, Phys. Rev. B 59, 14157 (1999). 

49 D. A. Huse, Phys. Rev. B 37, 2380 (1988). 

50 B. Efron and R. Tibshirani, An Introduction to the Boot- 
strap (Chapman & Hall/CRC, Boca Raton, 1993). 

51 J. D. Reger, J. A. Riera, and A. P. Young, J. Phys.: Con- 
dens. Matter 1, 1855 (1989). 

52 P. L'Ecuyer, Math. Comput. 65, 203 (1996). 

53 Although degeneracy is removed by disorder, we have to 
handle it if we want our algorithm to be valid in the undi- 
luted case. 

54 Subroutine DHSEQR already gives the product Y'Z in- 
stead of Z. 

The negligence of normalization doesn't change the conclu- 
sions we will arrive. Actually, this modes will be identified 
with the Goldstone modes of the diluted system, which, as 
in the clean limit, can have divergent amplitude. 

56 This procedure is a highly inefficient one, because only 
a small fraction (~ 6% for L — 14) of disordered lat- 
tice realization are accepted. Nevertheless, the amount of 
time spent finding clusters with N a — Nb is a small per- 
centage (~ 15% and ~ 6% for L — 14 in the Cholesky 
decomposition method and Bogoliubov- Valatin transforma- 
tion method, respectively) of the time consumed by the 
diagonalization subroutines. 

57 As the simulated lattices have N — 2 x L x L sites, and the 
dilution is achieved generating a random number r £ [0, 1] 
at each lattice site, we need to be careful with the period 
of the random number generator when averaging over 10 
disorder realizations. In this work we have used the max- 
imally equidistributed combined Tausworthe generator 
as implemented in the GNU Scientific Library. The period 
of this generator is 2 88 (~ 10 26 ). 

58 The Neel temperature determined from Eq. 1901 for 
isotropic Heisenberg systems is known to be overestimated. 
In order to correct for it a self consistent solution of this 
equation is used, where S is replaced by m av . This proce- 
dure, very simple to implement in the undiluted case, since 
there is an analytical expression available for the magnon 
spectrum, is much more difficult in our case. We postpone 
the discussion of this aspect to a latter publication. 



